Serum_Autoantibody_2021

library(car)
library(lsmeans)
library(calibrate)
library(dplyr)
#library(DEGreport)
library(DESeq2)
library(DEFormats)
library(edgeR)
library(emmeans)
library(ggpubr)
library(ggsci)
library(ggplot2)
library(gridExtra)
library(pheatmap)
library(reshape2)
library(RColorBrewer)
library(scales)
library(rstatix)
library(tidyr)
library(magrittr)
library(PCAtools)
library(tidyverse)
library(Biobase)
library(marray)
library(limma)
library(gplots)

0 Dataset

data <- read.csv("../../SerumAntibodies_2021/IgG_MCF_SH_179_raw.csv",  header=T, row.names = 1)
#Columns 16-47 are empty and are therefore removed
data <- as.data.frame((data))[,-c(16:47)]

#Creating metadata including Strain and age
Strain <- c(rep(c("BALBc", "NOD"), each=6), rep(c("NOD"), each=3))
Age <- c(rep(c("8wks", "28wks"), each=3), rep(c("8wks", "24wks", "33wks"), each=3))
Group <- c(rep(c("A","B","C","D","E"), each=3))
colData <- as.data.frame(cbind(c('B1_8', 'B2_8', 'B3_8', 'B1_28', 'B2_28', 'B3_28', 'N1_8', 'N2_8', 'N3_8', 'N1_24', 'N2_24', 'N3_24', 'N1_33', 'N2_33', 'N3_33'), Strain, Age, Group))
colnames(colData) <- c('Sample', "Strain", "Age", "Group")
rownames(colData) <- colData$Sample
colnames(data) <- colData$Sample

colData$Strain <- factor(colData$Strain)
colData$Strain <- relevel(colData$Strain, ref = "BALBc")
colData$Age <- factor(colData$Age, levels = c("8wks", "24wks", "28wks", "33wks"))
colData$Age <- relevel(colData$Age, ref='8wks')

1 Visualizing raw data

#Boxplot of raw Signal Distribution, not normalized

mydata<-as.matrix(log2(data+0.5))

boxplot(as.data.frame(mydata),main="Signal distribution, Not normalized",xlab="Samples",ylab="Log2 Intensity",cex=0.8)

######2 Clustering samples based on Pearson correlation ############

matcor<-cor(mydata)
range(as.vector(matcor))

## [1] 0.8503868 1.0000000

strain<-as.character(as.numeric(colData$Strain))
age<-as.character(as.numeric(colData$Age))

library(gplots)
heatmap.2(matcor,trace="none",col=heat.colors(40),ColSideColors=strain,
          RowSideColors=age, cexCol=1,cexRow=1,labCol="")

3 Data filtering and Normalization

mydata <- as.matrix((data))
dataN <- (mydata[1:95,]/mydata[96,] * 100) #df with all rows
dataN <- mydata
boxplot(as.data.frame(log2(dataN+0.5)),main="Signal distribution, Not normalized",xlab="Samples",ylab="Log2 Intensity",cex=0.8)

library(matrixStats) #for function rowmedians
IgG_SNR <- read.csv("../../SerumAntibodies_2021/IgG_MCF_SNR.csv",  header=T, row.names = 1)[,1:15]
IgG_SNR$average <- rowMeans(IgG_SNR)
IgG_SNR$med <- rowMedians(as.matrix(IgG_SNR))
mydata <- data[which(IgG_SNR$average>3),] #filtering rows with high noise (i.e., rows with Signal to Noise Ratio or SNR > 3)
mydata<-as.matrix(mydata)
dataN <- mydata
boxplot(as.data.frame(log2(dataN+0.5)),main="Signal distribution, Not normalized SNR>3",xlab="Samples",ylab="Log2 Intensity",cex=0.8)

dataN <- (mydata[1:69,]/mydata[70,]) * 100000
boxplot(as.data.frame(log10(dataN+0.5)),main="IgG normalization",col=strain)

library(limma) #for function plotMDS
plotMDS(log2(dataN+0.5))

mm <- model.matrix(~0 + Group)
y <- voom((dataN), mm, plot = T)

countData = as.data.frame(dataN)

df_dseq = melt(countData, variable.name = "Samples", value.name = "count") # reshape the matrix 
df_dseq = data.frame(df_dseq, Condition = substr(df_dseq$Samples,1,1))
mycolors <- colorRampPalette(brewer.pal(8, "Set1"))(13)

ggplot(df_dseq, aes(x = log10(count+0.5), colour = Samples, fill = Samples)) + 
   geom_density(alpha = 0, size = 0.8) + 
   facet_wrap(~Condition, ncol=3) #+

   #theme_minimal() + xlim(-1.5,6) +
   #scale_colour_manual(values=mycolors, name="") + guides(fill="none")

matcor<-cor(dataN)
range(as.vector(matcor))

## [1] 0.1740427 1.0000000

strain<-as.character(as.numeric(colData$Strain))
age<-as.character(as.numeric(colData$Age))
heatmap.2(matcor,trace="none",col=heat.colors(40),ColSideColors=strain,RowSideColors=age, cexCol=1,cexRow=1,labCol="")

4 DGE analysis

#mydata <- as.matrix((data+0.5))

conditions<- paste(colData$Strain[1:12], colData$Age[1:12],sep=".")

conditions <- factor(conditions, levels=unique(conditions))
design <- model.matrix(~0+ conditions)
colnames(design) <- levels(conditions)

#contrasts of the previous analysis. 
#cont.matrix<- makeContrasts(B8vsN8 = NOD.8wks - BALBc.8wks, B28vN24 = NOD.24wks - BALBc.28wks,
#   B28vN33 = NOD.33wks - BALBc.28wks,
#   N24vN8 = NOD.24wks - NOD.8wks,
#   N33vN8 = NOD.33wks - NOD.8wks,
#   N33vN24 = NOD.33wks - NOD.24wks,
#   levels = design)
#Introducing new contrasts by removing the 33 week old samples.
fit <- lmFit(dataN[,1:12], design)
cont.matrix<- makeContrasts(
B8vsN8 = NOD.8wks - BALBc.8wks,
B28vN24 = NOD.24wks - BALBc.28wks,
N24vN8 = NOD.24wks - NOD.8wks,
B28vB8 = BALBc.28wks - BALBc.8wks,
levels = design)

fit.cont<- contrasts.fit(fit, cont.matrix)
fit.cont<- eBayes(fit.cont)

decide <- matrix(c("fdr",0.05,
                   "fdr",0.1, "none",0.005, "none", 0.01),nrow=4,ncol=2,byr=T)
mysum <- as.list(1:nrow(decide))
mynum <- 0
maxmax <- 0
for (test in 1:nrow(decide)){
   results<-decideTests(fit.cont,
                        adjust.method=decide[test,1],
                        p=as.numeric(decide[test,2]))
   summary(results) -> mysum[[test]]
   mynum[test] <-length(which(apply(results,1,function(x)any(x,na.rm=T))))
   maxmax <- max(c(maxmax, as.vector(mysum[[test]][c(1,3),])))
}
par(mfrow=c(1,nrow(decide)))
for (test in 1:nrow(decide))
   {
   as.numeric(as.vector(mysum[[test]][3,]))->plotMe1
   as.numeric(as.vector(mysum[[test]][1,]))->plotMe2
   maxData = max(plotMe1)
   maxData2 = max(plotMe2)
   barplot(plotMe1,horiz=T,col="red",xlim=c(-maxmax,maxmax),
           main=paste("Gene Changes \np<",decide[test,2], ", " , decide[test,1],
                      " (" ,mynum[test] ,")",sep=""))->yy
   barplot(-plotMe2,horiz=T,col="green",add=T)->yy
   xx<-vector("integer",ncol(mysum[[test]]))
   text(xx,yy,colnames(mysum[[test]]))
   text((plotMe1+10)*0 + .9*maxData,yy+0.1,format(plotMe1,digits=3))
   text((-plotMe2-10)*0 - .9*maxData2,yy+0.1,format(plotMe2,digits=3))
}

results<-decideTests(fit.cont,adjust.method="none", p=0.01)
summary(results)

##        B8vsN8 B28vN24 N24vN8 B28vB8
## Down        0       0      0      0
## NotSig     68      65     67     69
## Up          1       4      2      0

write.fit(fit.cont,file="IgG_DGE_0423.csv",adjust="none",results=results, sep=',')
fitObj <- read.csv("IgG_DGE_0423.csv", header=T, row.names = 1)

toptable1 <- topTable(fit.cont, coef=2, adjust.method="none", p=0.01, sort.by = "P", n = Inf)
write.csv(toptable1,file="IgG_DGE_0423_II.csv")

myNames<-names(fitObj)
res.col<- which(regexpr("Res.",myNames)>0)
DElist<- which(apply(fitObj[,res.col],1,function(x)any(x,na.rm=T)))
length(DElist)

## [1] 5

fitObjDE <-fitObj[DElist,]

fitObjDE$Autoantigen <- rownames(fitObjDE)
fitObjDE

##                AveExpr Coef.B8vsN8 Coef.B28vN24 Coef.N24vN8 Coef.B28vB8
## Collagen VI 344353.360  306988.513    698156.44  466614.955  75447.0273
## Myosin       34663.933   68677.811     42702.87  -21566.851   4408.0887
## PCNA         25915.298   16328.886     79655.24   64814.316   1487.9658
## Sm            6144.959    6763.382     14148.64    7357.128    -28.1265
## ssDNA       753270.453  679724.733   1795054.96 1226257.802 110927.5744
##              t.B8vsN8 t.B28vN24  t.N24vN8     t.B28vB8 P.value.B8vsN8
## Collagen VI 2.6252712  5.970419  3.990347  0.645199725    0.028906317
## Myosin      4.6927709  2.917897 -1.473668  0.301205734    0.001322884
## PCNA        0.9255585  4.515041  3.673823  0.084341295    0.380216625
## Sm          1.7581652  3.677988  1.912512 -0.007311583    0.114578304
## ssDNA       1.4040943  3.708010  2.533057  0.229140955    0.195776818
##             P.value.B28vN24 P.value.N24vN8 P.value.B28vB8         F
## Collagen VI     0.000264043    0.003544981      0.5358308 17.760699
## Myosin          0.018153085    0.176632442      0.7704991 10.407852
## PCNA            0.001686093    0.005656377      0.9347390  9.434717
## Sm              0.005621207    0.090039827      0.9943351  6.144546
## ssDNA           0.005374459    0.033497880      0.8241801  6.511897
##                F.p.value Results.B8vsN8 Results.B28vN24 Results.N24vN8
## Collagen VI 0.0005205418              0               1              1
## Myosin      0.0032778110              1               0              0
## PCNA        0.0044953318              0               1              1
## Sm          0.0162032926              0               1              0
## ssDNA       0.0137576330              0               1              0
##             Results.B28vB8 Autoantigen
## Collagen VI              0 Collagen VI
## Myosin                   0      Myosin
## PCNA                     0        PCNA
## Sm                       0          Sm
## ssDNA                    0       ssDNA

library(devtools)
#install_github("dpgaile/AutoAntArrayExmpl")
library(AutoAntArrayExmpl)

library(devtools)
#install_github("dpgaile/AutoAntArrayExmpl")
library(AutoAntArrayExmpl)
library(NMF)
library(quantreg)
library(asbio)
library(fdrtool)
#library(discreteMTP)

#Sample Spot and Background Distributions

IgG_NSI <- read.csv("../../SerumAntibodies_2021/IgG_MCF_NSI.csv",  header=T, row.names = 1)[1:96,1:12]
IgG_SNR <- read.csv("../../SerumAntibodies_2021/IgG_MCF_SNR.csv",  header=T, row.names = 1)[1:96,1:12]
Strain <- c(rep(c("BALBc", "NOD"), each=6))
Age <- c(rep(c("8wks", "28wks"), each=3), rep(c("8wks", "24wks"), each=3))
Group <- c(rep(c("A","B","C","D"), each=3))
colData <- as.data.frame(cbind(c('B1_8', 'B2_8', 'B3_8', 'B1_28', 'B2_28', 'B3_28', 'N1_8', 'N2_8', 'N3_8', 'N1_24', 'N2_24', 'N3_24'), Strain, Age, Group))
colnames(colData) <- c('Sample', "Strain", "Age", "Group")
rownames(colData) <- colData$Sample
colnames(IgG_NSI) <- colData$Sample
colnames(IgG_SNR) <- colData$Sample

colData$Strain <- factor(colData$Strain)
colData$Strain <- relevel(colData$Strain, ref = "BALBc")
colData$Age <- factor(colData$Age)
colData$Age <- relevel(colData$Age, ref='8wks')

######Data Filtering #######

colData <- colData[1:12,] #remove the last group - M NOD 33 wks
IgG_SNR$average <- rowMeans(IgG_SNR)
IgG_SNR$med <- rowMedians(as.matrix(IgG_SNR))
IgG_NSI <- IgG_NSI[which(IgG_SNR$average>3),] #filtering rows with high noise (i.e., rows with SNR > 3)
IgG_SNR <- IgG_SNR[which(IgG_SNR$average>3),]
IgG_raw=list()
IgG_raw$NSI <- as.matrix(IgG_NSI)
IgG_raw$SNR <- as.matrix(IgG_SNR)[,1:12]
IgG_raw$SInfo <- colData

clrs=c(rep(pal_jco("default")(5), each=3))[1:12]
pchs=c(rep(3,3),rep(4,3), rep(5,3), rep(6,3), rep(7,3))[1:12]

#tiff("fig1_.tiff", units="in", width=7, height=7.4, res=300)
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
plot(as.vector(IgG_raw$NSI),as.vector(IgG_raw$SNR),type="n",ylab="IgG Background Intensity (SNR)",xlab="IgG Spot Intensity (NSI)" )
abline(0,1,lwd=3,col="steelblue")
for(j in 1:12) points(IgG_raw$NSI[,j],col=clrs[j],IgG_raw$SNR[,j], pch=j,cex=0.5)
plot(log(as.vector(IgG_raw$NSI)),
     log(as.vector(IgG_raw$SNR)),
     type="n",ylab="log IgG Background Intensity (SNR)",xlab="log IgG Spot Intensity (NSI)")
abline(0,1, lwd=3,col="steelblue")
for(j in 1:12) points(log(IgG_raw$NSI[,j]),log(IgG_raw$SNR[,j]),
                      pch=j,col=clrs[j],cex=0.5)
legend(6.5,0,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)

plot(density(log10(IgG_raw$NSI[,1]+ 0.5),type="l"), main="  log10(NSI +0.5)", sub=" . ")
for(j in 2:12) {
   dens<-density(log10(IgG_raw$NSI[,j]+0.5))
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(-0.8,0.6,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
plot(density(log10(IgG_raw$NSI[,j]+0.5)), main="log10(SNR +0.5) ", sub=" . ", type="l")
for(j in 1:12) {
   dens<-density(log10(IgG_raw$NSI[,j]+0.5))
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}

#dev.off()

#Centering and Scaling Differences Across Features

# Tukey Tri-mean for location
TriMnG=RowTriMeans(IgG_raw$NSI)
# Biweight midvariance for spread
bw.estG=unlist(apply(IgG_raw$NSI,1,r.bw))

#tiff("fig2.tiff", units="in", width=7, height=3.6, res=300)
par(mfrow=c(1,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,3.25,1.5,0.75))
plot(TriMnG,bw.estG,pch=c(rep(12)),
     col=c(rep("darkgreen",12)),
     xlab="Autoantigen Tukey TriMean",
     ylab="Autoantigen Biweight Midvariance")
legend("topleft",pch=c(12),col=c("darkgreen"),
       legend="IgG",bty="n")
plot(log(TriMnG),log(bw.estG),pch=c(rep(12)),
     col=c(rep("midnightblue",12)),
     xlab="Log Autoantigen Tukey TriMean",
     ylab="Log Autoantigen Biweight Midvariance")
legend("topleft",pch=c(12),col=c("midnightblue"),
       legend="IgG",bty="n") 

#dev.off()

#Heatmaps heatmaps using the uniform rank-its (i.e., rank values scaled to the unit interval) of the assays values across samples (i.e., the rank-its of the spot quantifications for each autoantigen across all 10 samples)

if(!file.exists("docs/Fheatmap3.png")){
   nmf.options(grid.patch=TRUE)
   RnkSmplW=matrix(nrow=70,ncol=12)
   RnkSmplW[,1]=IgG_raw$NSI[,1]
   for(j in 2:12){ RnkSmplW[,j]=IgG_raw$NSI[,j]}
   
   for(i in 1:70) {RnkSmplW[i,]=rank(RnkSmplW[i,])/(13)} #(row ranks)
   colnames(RnkSmplW)=as.character((colData$Sample))
   rownames(RnkSmplW) <- rownames(IgG_raw$NSI)
   
   sink("mess.txt")
   aheatmap(RnkSmplW,annCol=factor(colData$Strain),
            main="Rank-it (Across Samples) Raw Signal IgG",
            dist="euclidean",hclustfun="ward.D",
            filename="Fheatmap3.png")
   sink()
}
library(pheatmap)

#tiff("fig3_.tiff", units="in", width=6, height=12.5, res=300)   
pheatmap(RnkSmplW, annCol=factor(colData$Strain), main="Rank-it (Across Samples) Raw Signal IgG", dist="euclidean", hclust="ward")

#dev.off()

Normalization First Stage

eps=.10
#tiff("fig4.tiff", units="in", width=6, height=8, res=300)  
par(mfrow=c(4,3));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
TriMn=RowTriMeans(IgG_raw$NSI)
ResMat=(IgG_raw$NSI)-matrix(rep(TriMn,12),ncol=12)
SubMat=abs(ResMat)<quantile(as.vector(abs(ResMat)),prob=1-eps)
for(j in 1:12){
   plot(TriMn,ResMat[,j],
        xlab="Raw Signal Tri-Mean",ylab="Residual Raw Signal",
        main=paste("IgG Raw Signal ; Sample",j),
        cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
   y=ResMat[,j]
   x=TriMn
   fitj <- rq(y ~ x, tau = .5,subset=which(SubMat[,j])) #calculating the the quantile reggression fit, excluding the outliers
   abline(coef(fitj),col=clrs[j])
}

#dev.off() ### here we are plotting the autoantigen rowmeans as a function of their residuals.

enable to background correct

IgG_raw$NSI <- as.matrix(IgG_NSI)
IgG_raw$NSI=IgG_raw$NSI*IgG_raw$SNR
par(mfrow=c(1,1));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
plot(density(log10(IgG_raw$NSI[,1]+ 0.5)), main="  log10(NSI +0.5)", sub=" . ",type="l")
for(j in 2:12) {
   dens<-density(log10(IgG_raw$NSI[,j]+0.5))
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}

IgG_raw$X_norm_1=IgG_raw$NSI
eps=.10
IgG_raw$QregFits=array(NA,dim=c(12,70))
TriMn=RowTriMeans(IgG_raw$X_norm_1) ### rowmeans
ResMat=IgG_raw$X_norm_1-matrix(rep(TriMn,12),ncol=12)
SubMat=abs(ResMat)<quantile(as.vector(abs(ResMat)),prob=1-eps)

#tiff("fig5_.tiff", units="in", width=7, height=6.5, res=300)  
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
plot(rep(TriMn,12),as.vector(ResMat),
     xlim=quantile(TriMn,prob=c(0,.8)),
     ylim=quantile(as.vector(ResMat),prob=c(0.05,.95)),
     type="n",xlab="Raw Signal Tri-Mean",ylab="Raw Signal",
     main="IgG Raw Signal")
for(j in 1:12){
   points(TriMn,ResMat[,j],cex=0.5,pch=pchs[j],col=clrs[j])
   y=ResMat[,j]
   x=TriMn
   fitj <- rq(y ~ x, tau = .5,subset=which(SubMat[,j]))
   abline(coef(fitj),col=clrs[j])
   IgG_raw$QregFits[j,]=coef(fitj)
   prdct=coef(fitj)[1]+coef(fitj)[2]*TriMn
   IgG_raw$X_norm_1[,j]=IgG_raw$X_norm_1[,j]-prdct #model fit
}
AAplotDens(IgG_raw,SampleIndices=1:12,clrs=clrs,useMat="norm",
           xlims=NA,returnModes=T,main="Normalized IgG Spot Signal")

##  [1] 4256.431 4327.828 4185.033 4327.828 3970.842 3899.444 4185.033 2614.294
##  [9] 3613.856 3828.047 3114.075 3899.444

legend(2000000,0.0006,lty=1,lwd=3,col=clrs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)

plot(density((IgG_raw$QregFits[1,]),type="l"), main="(QregFits)", sub=" . ")
for(j in 2:12) {
   dens<-density((IgG_raw$QregFits[j,]))
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
plot(density(log(IgG_raw$X_norm_1[,1]+max(IgG_raw$X_norm_1[,1]))-min(log(IgG_raw$X_norm_1[,1]))), main="(X_norm_1) ", sub=" . ",type="l")
for(j in 2:12) {
   dens<-density(log(IgG_raw$X_norm_1[,j]+max(IgG_raw$X_norm_1[,1]))-min(log(IgG_raw$X_norm_1[,1])))
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(400000, 0.00005,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)

#dev.off()

#tiff("fig6.tiff", units="in", width=7, height=6, res=300)  
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
eps=0.10
q.eps=0.10
lambda=.10 # incremental change
nitr=95 # number of iterations (going upto 95, 46th iteration was found to be optimum)
sad=rep(0,nitr) # sum of absolute deviations

# re-init
IgG_raw$X_1step_norm_1=IgG_raw$X_norm_1 
IgG_raw$X_norm_1=IgG_raw$NSI

for(i in 1:nitr){
     # snoop and grab invariants
   X=IgG_raw$X_norm_1
   iWRS=function(i) t.test(X[i,1:6],X[i,7:12])$p.value
   qMat=matrix(rep(mapply(iWRS,1:70),12),ncol=12)
   qMat[which(is.na(qMat[,1])==T),] <- 0
   qcut=quantile(qMat[,1],prob=q.eps)
   
   TriMn=RowTriMeans(IgG_raw$X_norm_1)
   ResMat=IgG_raw$X_norm_1-matrix(rep(TriMn,12),ncol=12)
   SubMat=abs(ResMat)<quantile(as.vector(abs(ResMat)),prob=1-eps)
   SubMat[qMat<qcut]=FALSE
   for(j in 1:12){
      y=ResMat[,j]
      x=TriMn
      fitj <- rq(y ~ x, tau = .5,subset=which(SubMat[,j]))
      prdct=coef(fitj)[1]+coef(fitj)[2]*TriMn
      IgG_raw$X_norm_1[,j]=IgG_raw$X_norm_1[,j]-lambda*prdct #calculating actual values
      sad[i]=sad[i]+sum(abs(2*lambda*prdct))
   }
} #### the lowest sad or "sum of deviations" occurs at the 46th iteration.
#qcut at this i=46 is 0.0054

plot(1:nitr,sad,main="iterative adj IgG")
smoothScatter(IgG_raw$X_1step_norm_1,IgG_raw$X_norm_1, nbin = 300)
abline(0,1,col="steelblue")

plot(density((IgG_raw$X_1step_norm_1[,1]),type="l"), main="(X_1step_norm_1) ", sub=" . ")
for(j in 2:12) {
   dens<-density((IgG_raw$X_1step_norm_1[,j]))
   lines(dens, cex=1.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(400000, 0.00005,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)

plot(density((IgG_raw$X_norm_1[,1]),type="l"), main="(X_norm_1) ", sub=" . ")
for(j in 2:12) {
   dens<-density((IgG_raw$X_norm_1[,j]))
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}

legend(4000,0.0009,lty=1,lwd=2,col=clrs[c(1,4,7,10)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)

#dev.off()

Second Stage of Normalization

The normalized assay values are log transformed for purposes of variance stabilization and ease of visualization. Prior to the log transform, a constant is added to the normalized data such that all normalized assay values are greater than or equal to 1.

deltaG=min(as.vector(IgG_raw$X_norm_1))
IgG_raw$W_1=log(IgG_raw$X_norm_1-deltaG+1)
IgG_raw$W_2=log(IgG_raw$X_norm_1-deltaG+1)

#tiff("fig7.tiff", units="in", width=8, height=7, res=300)  
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
AAModesIgG=AAplotDens(IgG_raw,SampleIndices=1:12,clrs=clrs,useMat="resid",xlims=NA,returnModes=T,main="IgG Residual (Stage 1) log Spot Signal")
legend(-0.6,14,lty=1,lwd=2,col=clrs[c(1,4,7,10)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1) #did not plot

# get resids..
W_1=IgG_raw$W_1
RW=RowTriMeans(W_1)
# get residual matrices
RW_MAT=matrix(rep(RW,dim(W_1)[2]),ncol=dim(W_1)[2])
W_1=W_1-RW_MAT #value - trimean

W=W_1
bw.est=unlist(apply(W,2,r.bw)) ####bivariate midvariance is something like pearsons r2. variances across columns

# now, rescale
for(j in 1:12){
   W_1[,j]=W_1[,j]/sqrt(bw.est[j]) #sq. root of bivariance is an estimate of sd
}
# now, get new values..
IgG_raw$W_1=RW_MAT+W_1
IgG_raw$W_2=RW_MAT+W_1

AAplotDens(IgG_raw,SampleIndices=1:12,clrs=clrs,useMat='resid',
           xlims=NA,returnModes=T,main="IgG Residual (Stage 2) log Spot Signal")

##  [1] -0.003662881  0.201954756  0.300397177  0.346442825 -0.101311411
##  [6]  0.350412277 -0.405371468 -0.214043860 -0.314074062 -0.192608817
## [11] -0.298196252 -0.137036483

legend(2.5,0.55,lty=1,lwd=2,col=clrs,legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)

plot(density(log2(IgG_raw$W_1[,1]+65)-7.20), main="(logW_1) ", sub=" . ", type="l")
for(j in 2:12) {
   dens<-density(log2(IgG_raw$W_1[,j]+65)-7.20)
   lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}

AAplotByMean(IgG_raw,SampleIndices=1:12,clrs=clrs,pchs=pchs,MeanSubset=1:12,useMat="log",main="IgG log normalized")

AAplotByMean(IgG_raw,SampleIndices=1:12,clrs=clrs,pchs=pchs,MeanSubset=1:12,useMat="resid",main="IgG residuals log normalized")
#dev.off()